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Abstract 

The text of Laplace, Sur V application du calcul des probabilites a la philosophie naturelle, (Theorie 
Analytique des Probabilites. Troisieme Edition. Premier Supplement), 1820, is quoted in the context of 
the Gram-Schmidt algorithm. We provide an English translation of Laplace's manuscript (originally in 
French) and interpret the algorithms of Laplace in a contemporary context. The two algorithms given 
by Laplace computes the mean and the variance of two components of the solution of a linear statistical 
model. The first algorithm can be interpreted as reverse square-root-free modified Gram-Schmidt by row 
algorithm on the regression matrix. The second algorithm can be interpreted as the reverse square-root- 
free Cholesky algorithm. 

1 Introduction 

This translation work is inspired by the one of Pete Stewart [3] who translates from Latin to English the 
"Theoria Combinationis Observationum Erroribus Minimis Obnoxiae" of Gauss in the SIAM book "Theory 
of the Combination of Observations Least Subject to Errors, Part One, Part Two, Supplement." Stewart 
translates 101 original pages of Gauss, and he also provides an important contribution (28 pages) to place 
the work of Gauss in a historical framework. This manuscript is a more modest contribution. I translate 
thirteen pages and explain the relation of Laplace's algorithm with our contemporary algorithms. I would 
like to thanks Pete Stewart to have inspired me by his work. I also would like to thank Ake Bjorck for giving 
me my first version of Laplace's manuscript back in 2004 and Serge Gratton for useful comments on an 
early draft of the manuscript. 

The goal of Laplace is to compute the mass of Jupiter (or Saturn) from a system of normal equations 
provided by Bouvard and from this same system to compute the distribution of error in the solution assuming 
a normal distribution of the noise on the observations. The parameters of the noise distribution are not 
known. Laplace explains how to compute the standard deviation of two variables of a linear statistical model. 
His algorithm can be interpreted as performing the Cholesky factorization of the normal equations and then 
compute the two standard deviations from the Cholesky factor. A second method used by Laplace to justify 
the first is to perform a QR factorization of the regression matrix and compute the standard deviation from 
the R factor. Laplace was performing the QR factorization through the modified Gram-Schmidt algorithm. 

Laplace did not know what a factorization was, nor a matrix. I interpret his result through factorizations 
but certainly do not claim that Laplace invented all this. 

The first method which Laplace introduces consists in successively projecting the system of equations 
orthogonally to a column of the observation matrix. This action eliminates the associated variable from the 
updated system. Ultimately, Laplace eliminates all the variables but the one of interest in the linear least 
squares problem, which eliminates all the columns but one in the observation matrix. Laplace is indeed 
introducing the main idea behind the Gram-Schmidt algorithm (successive orthogonal projections.) Laplace 
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gives an example on a s-by-6 system. Once the observation matrix is reduced to a vector column, Laplace 
is able to relate the standard deviation of the variable of interest to the norm of this vector column and the 
norm of the residual vector. 

While Laplace could have stopped here and performed the modified Gram-Schmidt algorithm onto the 
original overdetermined system, he explains how to compute the norm of the projected column of observa- 
tions of interest directly from the normal equations. He observes that, if he performs a Cholesky factorization 
of the normal equations, the last coefficient computed will be equal to the norm of the last column orthogo- 
nally projected successively to the span of the remaining columns. In the mean time, Laplace observes that 
this method (Cholesky) provides a way to get the value of the solution from the normal equations. Laplace 
also generalizes this approach to more than one variable. 

Laplace has used the modified Gram-Schmidt algorithm as a tool to derive the Cholesky algorithm 
on the normal equations. Laplace did not interpret his results with orthogonality, in particular, he did not 
observe that, after orthogonal projection with respect to the last column, all the remaining projected columns 
were made orthogonal to that column. The orthogonal projections are interpreted as elimination conserving 
orthogonality with the residual. Laplace correctly explains and observes the property that all the remaining 
projected columns, after elimination/projection, are orthogonal to the residual of the least squares problem 
and that the residual vector is conserved. 

Laplace then uses his Cholesky algorithm to solve two 6-by-6 systems of normal equations given to 
him by the French astronom Bouvard to recompute the mass of Jupiter and Saturn, the originality of the 
work consists in assessing the reliability of these computations by estimating the standard deviation of the 
distribution of error in the solution. 

In Section [2] I set up the background for Laplace's work. This background is briefly recalled in Section 

1 of Laplace's manuscript. I chose not to translate this Section directly. It is failry hard to read indeed and I 
have preferred to explain it and refer to the equations in it. In Section [3] I provide a translation from French 
to English of Laplace's Sur V application du calcul des probabilites a la philosophie naturelle, (Theorie 
Analytique des Probabilites. Troisieme Edition. Premier Supplement), 1820. I have translatted Sections 2 
and 5 which represent pp.505-512 and pp. 516-520. 

This manuscript of Laplace is quoted in the book of Farebrother |f2j Chap.4] and the book of Bjorck [Q] 
p. 61]). In both books, the authors claim that Laplace is using the modified Gram-Schmidt algorithm. 

The text is available in French from the Bibliotheque Nationale de France website^] There is one typo (in 
the pages we are translatting). On p. 517, the seventh equation should read —668486", 70 = — 13208350z + 
413134432z' - 151992,0z'" - 34876,7 Z "'. 

We present some terminology used by Laplace. The overdetermined system of equations is named: les 
equations generates de condition des elements. The Linear Least Squares method is named: la methode la 
plus avantageuse. The poids (=weight) P of a normal distribution is related to the standard deviation a by 
P={2a 2 y l . 

2 Background 

The main goal of this manuscript is to provide a translation of Laplace's algorithmic contr ibution. However 
to put things into context, we start in this section with some background and notations. 

'See http://gallica.bnf . fr/ark:/12 14 8/bpt6k775 950. For Section 2, type 505 in the box Alter Page. For Section 5, 
type 5 16 in the box Aller Page. 
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2.1 Covariance matrix of the regression coefficients of a linear statistical model 

Laplace considers the classical linear statistical model 

Ax = b* +e, A G R sxn , b* e R s , 

where e is a vector of random errors, having normal distribution and we will denote a e the standard deviation 
of ||e||. In statistical language, the matrix A is referred to as the regression matrix and the unknown vector x 
is called the vector of regression coefficients. In our context, the matrix A is full rank. If e = (there is no 
error in the data), then we denote the solution of the consistent overdetermined linear system of equations 

as x*, 

Ax* = b*. 

Given a vector of s observations b, Laplace considers the linear estimate x given by the solution of the 
linear least squares method. We call e' the residual of the linear least squares solution 

e =Ax-b. 

In Laplace terms, (see, e.g., [p.501] last sentence), by the conditions de la methode la plus avantageuse, we 
have 

jyo e '(o = o, jyv(o = o, 

where pW is the element (i,l) of A, is the element (i,2) of A, . . ., and eW is the element i of e'. In other 
words, 

e' = Ax - b _L Span (A). 

Several other definitions for x, the linear least squares solution, are possible. We give two more equivalent 
definitions 

x is such that ||Ax — b|| = min ||Ay — b|| , or, equivalently, x is such that A T Ax = A r b. 

We define u = (u, u', . . .) the random vector which represents the error between the vector x* and the linear 
estimator x. From p. 501 to p.504, Laplace derives the formula of the joint distribution of the random 
variables u, u', ... 

We know that the joint distribution of the multivariate centered normal variable u is proportional to 



exp 



^C-'u), 



2 

where C is the covariance matrix of u. In our case we have 

therefore, the joint distribution of u is proportional to 

1 



exp(-^u r A r Au). (1) 



In practice one does not know <5 l h and we therefore rely on an unbiased estimate, for example 

— llb-Axll 2 . 



s — n 
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Now if we approximate s — n with s, we obtain that the random vector u follows a multivariate normal 
distribution with covariance matrix 

^||b-Ax|| 2 (A r A) _1 = ^lle'll 2 (A r A)^ . 

So that the joint distribution of the random variables u, u', ... is proportional to 



s 



exp {~2Mf uTATAu 

This result is given in term of the variable v = u/y/s by Laplace. Laplace states that the joint distribution of 
the random variables v, v', ... is proportional to (see first formula top of p.504) 

/ L ^(o v+ ^v + ...) 2 \ 

exp 



2££'(<)2 I • 

Note that 

«pf-^VA U W*p(- E(P< ' W ' V + '' )2 ' 



2||e'|| 2 y I 2£e'(02 / 

Therefore, Laplace's framework fits our standard linear statistical model framework. 



2.2 Laplace's algorithm to compute of the standard deviation of one variable of a linear 
statistical model 

2.2.1 Laplace and the modified Gram Schmidt algorithm 

The background is now half set. We have a linear statistical model to which we seek a regression vector x 
through the linear least squares method and we know that the covariance matrix of the regression vector is 
given by the matrix C. Laplace wants to compute only the first variable, z, of the regression vector, x. He 
also seeks the standard deviation, o z = o u , of this variable. 

From p.504 to p.505, Laplace explains that the modified Gram-Schmidt process applied to the matrix A 
enables him to find the standard deviation of the first variable. Laplace applies the modified Gram-Schmidt 
in a backward manner, that is, he projects the columns 1 to n — 1 orthogonally to the span of the column n 
and obtains the matrix Ai, then, working from the updated matrix Ai, he projects the columns 1 to n — 2 
orthogonally to the span of the column n — 1 , etc. 

The reverse modified Gram-Schmidt by row algorithm on the matrix A is formally given as follows. 

Following Laplace, we name the columns of A 

A = (p,q,r,...t,g,l). 
First step is to project column 1 to n — 1 of A 

Pi 

qi 



ti 

gi 



orthogonally to the span of its column n, so we define 
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This defines the s-by-(n — 1) matrix Ai = (pi,qi,ri, . . .ti,gi) . 

Second step is to project column 1 to n — 2 of A i orthogonally to the span of its column n — 1, so we 
define 




ti- 



Pi, 



qi, 



This defines the s-by-(n — 2) matrix A2 = (p2,q 2 ,i"2, . . .I2) . 

At the end of step n — 2, we have computed the s-by-2 matrix A„_2- The step n — 1 consists in projecting 
the first column of A„_2 orthogonally to the span of its second column 



Nowadays we are use to describe the modified Gram-Schmidt the other way around: project orthogo- 
nally to column 1 , then column 2, etc. In either case, we note that we need to order our variables correctly. 
With Laplace's method (reverse modified Gram-Schmidt), we will see that it is crucial to have the variable 
of interest ordered first. (And ordered last in the case of forward modified Gram-Schmidt.) 

Forward modified Gram-Schmidt generates a QR factorization of the matrix A, that is, we compute 

A = QR, where Q is s-by-n with orthonormal columns and R is n-by-n upper triangular. 

(Without loss of generality, we will impose the diagonal elements of R to be positive.) On the other hand, 
reverse modified Gram-Schmidt generates a QL factorization of the matrix A, that is, we compute 

A = QL, where Q is s-by-n with orthonormal columns and L is n-by-n lower triangular. 

(Without loss of generality, we will impose the diagonal elements of L to be positive.) 

We note that Laplace does not generate the matrix Q. Laplace generates the matrix T defined as 



If we normalize the columns of T, we will obtain Q. Laplace applies what we could call the reverse square- 
root-free modified Gram-Schmidt by row algorithm. If we define 




I = (p„_i,q„_2,r„_3,... 



t2,gl,l) ■ 



(2) 



D M = diag (||p„_i || 2 , ||q„- 2 || 2 , ||r„- 3 || 2 , • ■ ■ ||t 2 || 2 , ||gi || 2 , ||1|| 2 ) = T T T. 



(3) 



then we have 



T = QD 



1/2 



M 



So that we also have the factorization 




(4) 



The matrix I D M ' L J is lower triangular with ones on the diagonal. 
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QR factorization 



reverse square-root-free QR factorization 




2.2.2 A standard relation between the standard deviation of the last variable of a statistical model 
and the QR factorization of the regression matrix 

The standard deviation of the variable i of u is given by the square -root of the entry (/, i) of the covariance 
matrix C. If we are interested in the standard deviation a u of the first variable of u, u, we need to be able to 
compute the entry ( 1 , 1 ) of (A r A) . We outline below a standard way to compute this quantity. 
Once we have the QL factorization of A, we write 

A r A = I/Q/QL = L T L 

therefore 

a u = gentry (1,1) of C = a* gentry (1,1) of (A^A)- 1 = o b gentry (1,1) of L *L T . 
And, so using the fact that the matrix L is lower triangular, we have 

o u = ^-o b . (5) 
£i,i 

We can prove that 

^1,1 = ||Pn-l||j 

(where p„_i is the vector obtained at the last step of reverse modified Gram-Schmidt algorithm), so that we 
obtain that the marginal probability density function of the first variable z is 

CXP (~2c?" P ' 1 ~ 1 " 2m2 ) ' ^ 

and if we use the fact that - s ||e'|| 2 can be used as an approximation of an unbiased estimate of o|, we obtain 
that the marginal probability density function of the first variable z is 

1 s " '|2„2 



6XP V 2jeT U 

This formula is assessed by Laplace on top p.505. We read: 
"This exponential becomes 

exp (— Pu ) , 
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where 

„ (i)2 

w error of the random variable z, P is what I called the poids ( weight) of this value." 

The poids is related to the standard deviation a with P = (2a 2 ) 1 . The term poids was chosen by 
Laplace for the following reason (see p.499): 

'7c? probability decroit avec rapidite quand il [le poids] augmente, en sorte que le 
resultat obtenue pese, sije puis ainsi dire, vers la verite, d'autant que ce module est 
plus grand." 

which gives 

"the probability quickly decreases with it [le poids] increases, so that the result 
weights, if I can says so, towards the truth as much as this modulus is larger." 

Other reasons are given in the same paragraph. 



2.2.3 Laplace's derivation of the standard deviation of the last variable from the QL factorization of 
the regression matrix 

The overall strategy of Laplace to compute o u is well-known nowadays. How did Laplace derive it in the 
first place? Starting from the fact that the joint density function of u, u' , . . ., u^ n > is proportional to 

exp f — -^||Am|| 2 

(see Equation([T]), Laplace is interested in computed a function proportional to the marginal probability 
density function of the first variable u, G„, that is Laplace wants to compute a function of the variable u 
proportional to 

/ / .../ exp d|Au|| 2 )du'du" ...du [n) 

Ju'=-°°Ju"=-°° J «(»)=-« V a b / 

Laplace proposes to proceed by steps. First we will seek a function proportional to the joint density 
function of u, u' , . . ., i/" -1 ); then we will seek a function proportional to the joint density function of u, u', 
. . ., u( n ~ 2 \ etc. we will eventually end up with a function proportional to the marginal probability density 
function of the first variable u. 

To perform the first step, we therefore need to compute a function of the variables 
proportional to 



[ + °° expf - \\\ Anil 2 ] du^ n \ 
Jd«)=-o° V al / 



J b 

Laplace observes that, (Pythagorean theorem), 



Au|| 2 = II 1/-^ I Au|| z + ||^Au 



ll r \ . ,„ „ 11' 



hi 2 / imp 
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and so 



/ " \ 



= HA, 



/ u \ 



( u \ 



rl ? A 



/ " \ 



l|Ai 



| 2 + IIHI 2 



1 



HI 



l r (p,q,...,t,g) 



V 



The joint density function of u, u' , . . ., is therefore proportional to 

/ / 



/ / u \ \ 



exp 



1 



|Ai 



V 



•exp 



1 



111 



V 



/ u \ \ 2 \ 



1 



HI 



l r (p,q,...,t,g) 



V 



v / / y 



(This latter equation corresponds to the seventh equation on p.504.) 

As previously explained, the first step of Laplace's derivation consists in integrating this last term for 
u ( n ) ranging from — oo to +oo in order to obtain a function proportional to the joint density function of u, u' , 
. . ., u,( n ~ x \ So let us do this. We write 



r+°° 
•/«(«)=- 



exp 



/ u \ 



|Ai 



\ 



( I 



• exp 



" W + Jpl T (p,q,---,t,g) 



exp 



/ u \ 



1 



|Ai 



V 



The second term is of the form 



V v 



(»)=- 



exp 



" ( ' l) + ¥p lT(p ' q '---' t,8) 



V v 



/ u \\ 2 \ 



v « (n - 1} J J J 
( " \ \ 2 ^ 



v » {n - x) ) J ) 



du {n) 



du {n K 



exp (-n (u^ +g(u,u', . . . , M ("- 1 ))) 2 ^ du (n \ 

We note that this term is independent of the variables ^ K Therefore we can remove this term 

from the previous equation and conclude that the joint density function of u, u' , . . ., j/" -1 ) is proportional to 



/ 



/ u \ 



\ 



exp 



|Ai 



(7) 
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Continuing the process, we end up with the fact that the marginal probability density function of the first 
variable u is proportional to 




We recover Equation Q also given on top p. 505. 

From this equation, Laplace deduce that, to compute the standard deviation of u, he needs to compute 
||Pb-i II 2 - While it is clear that he (or Bouvard) can work on A and perform the modified Gram Schmidt 
algorithm, Laplace finds it easier to work on the normal equations. Quoting Laplace: 

"Mais il est plus simple d'appliquer le procede dont nous venons defaire usage awe 
equations finales qui determinent les elements, pour les reduire a une seule, ce qui 
donne une methode facile de resoudre ces equations :" 



which means 



"But it is easier to apply the method we have just used to the final equations which 
define the variables, in order to reduce it to a single, which gives a convenient way of 
solving these equations." 

Therefore the next question that needs to be answered is: how can we compute ||p«-i|| 2 from A T A 
without accessing A? This question is the matter of Section 2 of Laplace's treatise from p.505 to p. 5 12. An 
example of application of the technique is proposed in the Section 5 of the same manuscript from p.516 to 
p. 520. We provide a translation of these two parts in the next section. 

If we consider the QL factorization of the matrix A given by 




Ci, the covariance matrix of joint normal distribution of the variables u, u\,. . ., u n -\, is 

Ci=oi(A[Ai) _1 =of (L[Li)~ ! . 

We can derive this relation from Laplace's analysis for example. Another way to derive, this result is to 
remember that the covariance matrix Ci of the joint normal distribution of the variable u, u\,. . ., m„_i is the 
(n — l)-by-(n — 1) block of the covariance matrix C of the joint normal distribution of the variable u, u\,. . ., 
So if we write 

al (A r A) 



we see that the the (n — l)-by-(« — 1) block is a 2 , (L[Li) 1 = Ci. This is two ways to explain a standard 
result. 



T T 



z 

a 



LfLi +zz t 
az T 



Li 

z 1 ' 

za 

a 2 



r (L[L,; 



(MLi 



a 



-l 



LfLi) 



-l 



z r (L[Li) 2+1 
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3 Translation 

We now present a translation of Laplace's text. We proceed by couple of pages. First page gives the French 
version. Second page gives the translatted version. We recall that the notation S stands for £. 
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2. Reprenons l'equation generale de condition, et, pour plus de simplicite, bornons-la aux six elements 
z, z' , z", z'" , z' v , z v ; elle devient alors 

e« = p (0 z + g («V + r^z" + ? ( 'V" + + X (, V - a w . (1) 

En la multipliant par A,W et reunissant tous les produits semblables, on aura 

le signe integral S s'etendant a toutes les valeurs de i, depuis i = jusqu'a i = s — 1, s etant le nombre des 
observations employees. Par les conditions de la methode la plus avantageuse, on a Sk®£® = 0; l'equation 
precedente donnera done 

v _ iv SX^ m Sti?*t® „ sW f ) , Sk®q® Sl®pM SI® a® 

z ~~ z 1&M^~ Z IP) 2 " -2 ~Wj^~ z sW 2 ~ z sx(02 + sW 2 ' 

Si Ton substitue cette valeur dans l'equation ([I]) et si Ton fait 



on aura 





= yo: 




4° 


= f(0 




•P 


= r« 




J? 


= 




(0 

Pi 


= P { ' 






= a ( 





5^(0 r (0 

Sk®p® 
SA.02 ' 

S3LC02 ' 



e (0 = p (0 z + 9 fV + r (0 z // + ^ + Y (o z n. _ a (0 . 



(2) 



par ce moyen, 1' element z v a disparu des equations de condition que represente l'equation Q. En multipliant 
cette equation par y^ et reunissant tous les produits semblables, en observant ensuite que Ton a 

SyfV = 



en vertu des equations 

= Sk®e®, = S^>E® 
que donnent les conditions de la methode la plus avantageuse, on aura 



d'oii Ton tire 



,v ^tf'M ^yW ,Syfqf SyUp? , Syf'M 

Z = — Z 77T7, Z TTTT Z 7777, Z 



(02 



Sy\ 



,(02 



o (02 c (02 ' 

Syi ' Sy\ 1 
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2. We consider again the overdetermined system of equations, and, for the sake of simplicity, we restrict 
it to the six elements z, z', z", z'", z' v , z v ; it then becomes 

e (0 = p (0 z + ? («V + A'k" + + fV" + X®z v - a w . (1) 

Multiplying by A,W and grouping all similar products, we have 

Sk®e® = zSk®p® +lSk®q® + . . .-Sk®a®, 

the integral sign S ranging for all the values of i, from i = to i = s— 1, s being the number of observations. 
By the conditions of la methode la plus avantageuse, we have SX^E^ = 0; the former equation consequently 
gives 

v _ fv SX,('V0 ,„ SW) „ SW0 , S*/'V''> ^WgjO gXWrjO 
2 ~~ Z ~Sp)^~ Z 1P)^" Z ~Sp)2~~ Z SX(02 ~ Z SM02 + 5X('")2 ' 
If we replace this value in Equation ([T]) and if we perform 



Sk® 2 ' 

W _ r (0_ X (0^ ( ' V(0 



5X(«')2 

Pi = P 



(0 _ „(0_ X (0^ ( °9 (0 
(0 _ 

5X(02 ' 



SX(02 



we have 



e« = p W z + gtij + r f)/ + ,[<V" + ^ _ a (0 . (2) 

by this technique, the element z v has disappeared from the system of equations represented by Equation Q. 
Multiplying this equation by grouping all similar products, and observing that we have 

SYi°e (0 = 

from the equations 

= Sk®e®, = 5^e« 
given by the conditions of la methode la plus avantageuse, we have 

= zSyfpf +z'siSf +J'si?rf +z"'Sy { ; ) t[ i) +z h 'Sy { ? 2 -S$a®; 



from which we draw 



„,sjhl_ „s£rf_ ,S^ql_ Stfp® , SyW 
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Si Ton substitue cette valeur dans l'equation ([2]) et si Ton fait 

r 2 - h h _ (i) 2 



2 " 1 Yl 5 Y p ' 

c (0 (0 

p (0 (0 

Pi = Pi Yi — ^> 

s Y r 

u 2 — u t ^ ^2 



5y' 



on aura 



e (0 = p« z + ? « z / + r (0 z // + , (0^/ _ a f . (3) 
En continuant ainsi, on parviendra a une equation de la forme 

n®=pfz-a$- (4) 

II resulte du n°20 du Livre II que, si la valeur de z est determinee par cette equation et que u soit l'erreur de 
cette valeur, la probabilite de cette erreur est 



sSp 



25£ , (') 2 7I 



(02 ,s P f 

5 p 25e'W2 i 



,2 



Se'^ 2 etant la somme des carres des restes des equations de condition, lorsqu'on y a substitue les elements 

sSp m 

determines par la methode la plus avantageuse. Le poids P de cette erreur est done egal a 2 ^J {i)2 • 

II s'agit maintenant de determiner Sp^ 2 . Pour cela, on multipliera respectivement chacune des equations 
de condition representees par l'equation ([T]), d'abord par le coefficient du premier element, et Ton prendra 
la somme de ces produits; ensuite par le coefficient du second element, et Ton prendra la somme de ces 
produits, et ainsi du reste. On aura, en observant que par les conditions de la methode la plus avantageuse 
5p(0g(0 = 0, St/OgW = 0, . . ., les six equations suivantes : 



pa = 


p (2) 


z 


+ 


pq 


z 1 


+ 


pr 


z" 


+ 




z'" 


+ 


py 


z iv 


+ 


pX 


z v 


qa = 


pq 


z 


+ 


q® 


z 1 


+ 


qf 


z" 


+ 


qi 


z'" 


+ 




z iv 


+ 


qX 


z v 


7a = 


rp 


z 


+ 


fq 


z' 


+ 


r (2) 


z" 


+ 


ft 


z'" 


+ 


rf 


z iv 


+ 


7x 


z v 


ta = 


tp 


z 


+ 


tq 


z' 


+ 


Tr 


z" 


+ 


f(2) 


z'" 


+ 


ty 


z iv 


+ 


tX 


z v 


ya = 


JP_ 


z 


+ 


W 


z' 


+ 


yr 


z" 


+ 




z'" 


+ 


y(2) 


z iv 


+ 


ix 


z v 


Xa = 


Xp 


z 


+ 


Xq 


z' 


+ 


Xr 


z" 


+ 


fa 


z'" 


+ 


Xy 


z iv 


+ 


xw 


z v 



(A) 



oil Ton doit observer que nous supposons 

pW=Sp®\ pq = Sp%®, q®=SqW 2 , qr = Sq {i) r^, 



13 



If we replace this value in Equation (|2]) and if we perform 

h ~ r i 'l _ d)2 



JO _ JO JQ jW 

c (0 (0 

^2 - 9i ri „ (,) 2 ■ 



5y' 



o (0 (0 



o (0 (0 

a (0 _ a w_ y «w 
sir 



we have 



e (0 = p f z + + r « z " + t f{" - af . (3) 
Continuing in a similar manner, we end up with an equation of the form 

S= p f z -af. (4) 

From n°20 of Livre II, we know that, if the value of z is determined by this equation and if u is the error of 
the value, the probability of this error will be 



sS P f 

h e 2Se'W2 



where Se'^ 2 is the sum of the squares of the residuals of the equations of condition, after we replaced the 

sSp ( ' )2 

elements determined by la methode la plus avantageuse. he poids P of this error is then equal to 2S J^ 2 ■ 

Our next task is to determine Sp^ 2 . For this, we multiply each of these equations represented by 
Equation ([I]), first by the coefficient of the first element, and we take the sum of these products; then by the 
coefficient of the second element, and we take the sum of these products, and so on for the remaining. We 
have, by observing that the conditions of la methode la plus avantageuse S/>We(0 = q, SgWgW = 0, . . ., the 
six following equations: 



pa = 


p (2) 


z 


+ 


pq 


z' 


+ 


pr 


z" 


+ 


pt 


z'" 


+ 


py 


z iv 


+ 


pk 


z v 


qa = 


pq 


z 


+ 


qV) 


z' 


+ 


qf 


z" 


+ 


qt 


z'" 


+ 




z iv 


+ 


qk 


z v 


7a = 


Yp 


z 


+ 


fq 


z' 


+ 


r (2) 


z" 


+ 


ri 


z'" 


+ 


ry 


z iv 


+ 


7k 


z v 


ta = 


tp 


z 


+ 


tq 


z' 


+ 


tf 


z" 


+ 


t^ 


z'" 


+ 


ty 


z iv 


+ 


7k 


z v 


ya = 


JP_ 


z 


+ 


yq 


z' 


+ 


yr 


z" 


+ 


yt 


z'" 


+ 


y(2) 


z iv 


+ 


yk 


z v 


ka = 


kp 


z 


+ 


kq 


z' 


+ 


h- 


z" 


+ 


kt 


z'" 


+ 


ky 


z iv 


+ 


k& 


z v 



(A) 



where we have defined 

pW=Sp® 2 , pq = S P {i) q { % q®=Sq® 2 , qf = Sq {i) 7% 
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Si Ton multiplie pareillement les equations de condition representees par 1' equation Q respectivement 
par les coefficients de z et que Ton ajoute ces produits, ensuite par les coefficients de z' en ajoutant encore ces 



produits, et ainsi de suite, on aura le systeme suivant d'equations, en observant que S/?[ eW = 0, Sq^'B^ = 0, 
. . ., par les conditions de la methode la plus avantageuse, 



PlCCi 

qla[ 



yTaT 



pm 

Pin 

pih 
Mi 



+ 
+ 
+ 
+ 
+ 



pm 

(2) 

qln 
qlh 



+ 
+ 
+ 
+ 
+ 



pm z 



q\r\ 

(2) 

nh 
nyi 



+ 
+ 
+ 
+ 
+ 



pih 
qlh 
nh 

t {2) 
'i 



+ 
+ 
+ 
+ 
+ 



nTi 
hYi 

(2) 

Ti 



(B) 



ou Ton doit observer que 



pm — *>P\ ?i ) P\ — bp\ > 



En substituant, au lieu de p^\ q^ , leurs valeurs precedentes, on a 



ou 



on a pareillement 



piqi = Sp { >q 



(0„(0 



piqi =pq 



5P)2 



X(2) ' 



(2) 
Pi 

(2) 

Pin 



,(2) 
(2) 



9 
7?r 



Xp 
Xq~ 
Xp Xr 



Ainsi les coefficients du systeme des equations (|B]) se deduisent facilement des coefficients du systeme des 
equations ( [A] ). 

Les equations de condition representees par 1' equation Q donneront semblablement le systeme suivant 
d'equations 



P2M2 
qi^2 



t 2 a 2 



et Ton a 



(2) 
Pi 

Plqi 
pin 



Pih 



+ 
+ 
+ 
+ 



Piqi 

(2) 

q\> 
qlfi 
qlfi 



+ 
+ 
+ 
+ 



pm z 



(2) 

Pi 

plql 



Pi 



(2) 



pm 



qm 

r (2) 
^2 

YiPT 

(21 7 

Yj 

Y1P1 giVi 



Yi 



/2) 



Y1P1 Yi«i 

(2) 

Yi 



+ 
+ 
+ 
+ 



Pih 
qlh 
rih 

,(2) 



(C) 
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If we multiply similarly the equations represented by Equation Q respectively by the coefficients of z 
and we add these products, then by the coefficient of z' adding again these products, and so on, we have 



the following system of equations, by noting that Sp[ 
methode la plus avantageuse, 



0, Sq\ eW =0, . . ., from the conditions of la 



<7iai 
rial 



flOti 

YiaT 



(2) 
Pi 

piqi 
pm 

pih 

Mi 



+ 
+ 
+ 
+ 
+ 



pm 

(2) 
9l 

qih 



+ 
+ 
+ 
+ 
+ 



pm 
qWi 

(2) 

nh 



+ 
+ 
+ 
+ 
+ 



Piti 

qih 

t {2) 



+ 
+ 
+ 
+ 
+ 



nyi 
iiyi 

(2) 

Yi ; 



(B) 



where we have defined 



Substituting py ', qf\ 



or 



we have similarly 



Spfq^ 



(2) 
Pi 



Sp 



(02 



Pl9l 

with their previous values, we have 



plql = SpWq^ 



pxqx = pq- 



(2) 
Pi 

(2) 

pin 



p(2) 

7?r 



5X(')2 

T(2T ; 



Xp 2 
Xq 

Ml 

Xp Xr 



PlOCi 



^.p fan 

P5 ' 



Doing so, the coefficients of the system of equations (|B]) are easily computable from the coefficients of the 
system of equations ( |A| ). 

The equations represented by Equation ([3]) similarly give the following system of equations 



Pi^i 
qi^i 

Wi 



(2) 
Pi 

Piqi 
Pin 
Pih. 



+ 
+ 
+ 
+ 



Piqi 

(2) 
92 

^2 
^2?2 



+ 
+ 
+ 
+ 



W2 Z 



r 2 



z 



+ 
+ 
+ 
+ 



p 2 t 2 

q~lh 



(C) 



and we have 



(2) 
Pi 

Piqi 



p 2 a 2 



(2) 
Pi 

piqi 



PlO.1 



y\pr 

(2) ) 

Y1P1 giYi 



Yi 



Y1P1 Yi«i 

(2) 

Yi 
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On a pareillement le systeme d'equations 



P303 

^303 
^303 



P3 



(2) 



z 


+ 


P3q3 




+ 


z 


+ 


(2) 
?3 


z' 


+ 


z. 


+ 






+ 



+ ^3^3 

(2) 



(D) 



en faisant 



(2) 
Pi 

P3OC3 



(2) 
/>2 



P2<2 
,(2) > 

P2t2 qth 



f (2) 

t 2 p2 ?2«2 

re) • 

r 2 



on aura encore 



en faisant 



/}4«4 
^404 



(2) 

Pa 

PA^4 



Pa 



(2) , , 

Z + p\q\ Z , 



/?4<?4 Z + 



(2) 
P~3^3 



(2) 
<7 4 



P3''3 

(2) ) 

P3^3 43 r3 

(2) ; 

g _ 

P3>"3 «3Q 

,(2) • 



(E) 



Enfin on aura 
en faisant 



P5CX5 = pf ] z, (F) 



(2)_ (2) j?4ff4 _ P4^4 g4CC4 , 
P5 —Pa (2)~> P5&5— P4M4 W) • 



pf^ est la valeur Sp^ 2 , et le poids P sera 



<7 4 #4 



(2) 
^5 



25e'( ! ')2 " 

On voit par la suite des valeurs p( 2 \ p[ , p£\ ■ ■ ■ qu'elles vont en diminuant sans cesse, et qu'ainsi, pour le 
meme nombre d' observations, le poids P diminue quand le nombre des elements augmente. 

Si Ton considere la suite des equations qui determinent 775015, on voit que cette fonction, developpee 
suivant les coefficients du systeme des equations ([A]), est de la forme 

p~a + Mqa + Nra + . . . , 

le coefficient de p~a etant l'unite. II suit de la que si Ton resout les equations (|A]), en y laissant p~a, qa, 
fa, . . . comme indeterminees, -4^ sera, en vertu de l'equation (Fi, le coefficient de p~a dans l'expression 

de z.- Pareillement, -L sera le coefficient de qa dans l'expression de z!\ ~m sera le coefficient de roc dans 

l'expression de z"; et ainsi de suite du reste; ce qui donne un moyen de simple d'obtenir py , qy, . . .; mais 
il est plus simple encore de les determiner ainsi. 
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We similarly have the system of equations 



P3OC3 
^303 
^303 



P3 



(2) 



P3q3 

Wf3 



+ 
+ 
+ 



P3q3 

(2) 
03^3 



+ 
+ 
+ 



P3^3 Z 



q 3 r 3 

J 2 ) 



Z 



(D) 



by doing 



we also have 



by doing 



(2) 

P~3~q3 
P3U3 



(2) 
Pi 

P2U2 



,m > 

h 

f (2) ' 
12 

,P) ! 
r 2 



P4OC4 
q 4 a 4 



(2) 



(2) 
P\ 

P~4~q~4 



+ 
+ 



p 4 q 4 

(2) 

q\ 



z 1 , 
z!, 



(2) 



P3TT 
(2) ) 
r 3 

Pin gin 

(2) J 
r i 

«3^3 

(2) 



(E) 



Finally we have 
by doing 



(2) 

p 5 a 5 =p y 5 'z, 



(2) (2) P4q4 

q\ 



pf' is the value Sp^ 2 , and le poids P is 



p 5 a 5 = /J4OC4 



J 2 ) 



p 4 q 4 q4^4 



(2) 
04 



^5 
25£'(')2 ' 



(F) 



We see from the sequence of values // 2 \ , p^, • • • that they always go diminishing, and so, for the same 
number of observations, le poids P decreases when the number of elements increases. 

If we consider the sequence of equations which determine p$ 0C5, we see that this function, developed 
according to the coefficients of the system of equations dAl), is of the form 



pa+Mqa+Nra + . 



the coefficient of pa being the unity. It follows from there that if we solve the equations (|A]), by leaving pa, 
qa, fa, ... as unknowns, -L is, due to Equation (Fi, the coefficient of pa in the expression of z. Similarly, 



Ph 

L is the coefficient of qa in the expression of z'\ -m is the coefficient of fa in the expression of z"; and so 



(2) (2) 

on for the others; this gives a simple mean to obtain py , qy , . . .; but it is even simpler to compute them as 
follows. 
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D'abord 1' equation |f| donne la valeur de et de z. Si dans le systeme des equations |eJ) on elimine 
z au lieu de z', on aura une seule equation en z', de la forme 



(2),/. 



^ 5 a 5 = q" 5 'Z 

en faisant 



(2)_ (2) P4<?4 _ Pm P4^4 

<?5 -?4 f2T> ^as-^o^ T^j — • 

Pa P\ 

Si dans le systeme des equations (|d]) on elimine z au lieu de z", pour ne conserver a la fin du calcul que 
z" , on aura r$ en changeant dans la suite des equations qui, a partir de ce systeme, determinent la 
lettre p dans la lettre r, et reciproquement. On aura ainsi 

r (2) _ (2) _ pjr? 
4 _ r 3 (2) ) 

P3 

p 3 q 3 p 3 r 3 



r A q A = r 3 <73 (2) 

„(2) _ _(2) _ m? 

H\ — </3 p (2) ; 

r (2) _ (2) _ gsk 



Pour avoir t!p, on partira du systeme |c]), en changeant, dans la suite des valeurs de pf^, p^q-i, 
qYn,, . . ., la lettre p dans la lettre t, et reciproquement. 

(2) J 1 

On aura pareillement la valeur de y 5 , en partant du systeme des equations ( B 1 et changeant dans la suite 

(2) (2) 

des valeurs de p\ , p\ , . . ., la lettre p dans la lettre y, et reciproquement. 

(2) (2) (2) 

Enfin, on aura la valeur de X 5 en changeant, dans la suite des valeurs de p\ , p 2 , ■ ■ ■, la lettre p dans 

la lettre X, et reciproquement. 
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Firstly Equation |fJ) gives the value of and of z. If in the system of equations |eJ) we eliminate 
instead of z', we have a single equation in z', of the form 



(2),/. 



q 5 a 5 = q" 5 'Z 

by doing 



(2)_ (2) 7?4<?4 _ P4q4 P4^4 
<?5 —1a (2T' <?5a5 — <?4«4 7j\ • 

P\ Pa 

If in the system of equations {Dl we eliminate z. instead of z", in order to only keep at the end of the 
(2) — 

computation z", we have r 5 by changing in the sequence of equations which, starting from this system, 

(2) 

determine p 5 , the letter p by the letter r, and reciprocally. We then have 

(2) _ (2) _ 

r 4 — r 3 (2) > 

P3 

p 3 q 3 p 3 r 3 



r 4 q 4 = r 3 q 3 (2) 

J 2 ) J 2 ) _ ml 

H4 — H3 (2) > 

Pi 

r (2) _(2) _ P4g4 2 

r 5 — r 4 (2) > 

«4 



In order to have ^ , we start from the system (O, by changing, in the sequence of values of p 3 , P3qj, . . ., 
r 3 , q3>"3, ■ ■ ■> the letter /? by the letter t, and reciprocally. 

We similarly have the value of Y5 , starting from the system of equations (fei and changing in the 

(2) (2) 

sequence of values of p 2 , p 3 , . . ., the letter p by the letter y, and reciprocally. 

(2) (2) (2) 

Finally, we have the value of }i 5 by changing, in the sequence of values of p\ , p 2 , . . ., the letter p by 

the letter X, and reciprocally. 
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5. Appliquons maintenant cette methode a un exemple. Pour cela, j'ai profite de 1' immense travail que 
Bouvard vient de terminer sur les mouvements de Jupiter et de Saturne, dont il a construit des Tables tres 
precises. II a fait usage de toutes les oppositions observees par Bradley et par les astronomes qui l'ont suivi : 
il les a discutees de nouveau et avec le plus gran soin, ce qui lui a donne 126 equations de condition pour 
le mouvement de Jupiter en longitude et 129 equations pour le mouvement de Saturne. Dans ces dernieres 
equations, Bouvard a fait entrer la masse d'Uranus comme indeterminee. Voici les equations finales qu'il a 
conclues par la methode la plus avantageuse : 

7212", 600 = 79593 8z - 12729398z' + 6788, 2z" - 1959, 0z'" + 696, 13z ! ' v + 2602z v , 

-738297", 800 = - 12729398z + 424865729z' - 153 106, 5z" - 39749, lz'" - 5459z !V + 5722z v , 

237" ,782 = 6788 , 2z - 1 53 106, 5z' + 7 1 , 8720z" - 3 , 2252z"' + 1 , 2484z' v + 1,337 lz v , 

-40", 335 = - 1959, 0z- 39749, lz' -3,2252z" + 57, 1911z"' + 3,6213z ,v + 1, 1128z v , 

- 343" , 455 = 696 , 1 3z - 5459z' + 1 , 2484z" + 3 , 62 1 3z'" + 21, 543z' v + 46 , 3 1 0z v , 

- 1 002" , 900 = 2602z + 5722z' + 1 , 337 lz" + 1 , 1 1 28z'" + 46 , 3 1 0z' v + 1 29z v . 

Dans ces equations, la masse d'Uranus est supposee j^^; la masse de Jupiter est supposee z" 
est le produit de l'equation du centre par la correction du perihelie employe d'abord par Bouvard; z'" est la 
correction de l'equation du centre; z lv est la correction seculaire du moyen mouvement; z v est la correction 
de l'epoque de la longitude au commencement de 1750. La seconde du degre decimal est prise pour unite. 

Au moyen des equations precedentes renfermees dans le systeme ([A]), j'ai conclu les suivantes, ren- 
fermees dans le systeme (|B]) : 

27441",68 = 743454z- 12844814z' + 6761, 23z"- 1981, 45z w - 237, 97z ,v , 

-693812",58 = - 128448 14z + 4246 11920z' - 153165, 81z" - 39798, 46z'" -7513, 15z ! \ 

248", 1772 = 6761 , 23z - 153165, 81z' + 71 , 858 lz" - 3, 2367z'" + 0, 7684z'\ 

-31", 6836 = -1981, 45z - 39798, 46z' - 3 , 2367z" + 57,181 5z'" + 3 , 22 1 8z ,v , 

16", 5783 = -237,97z-7513,15z / + 0,7684z" + 3,2218z , " + 4,9181z ,v . 

De ces equations, j'ai tire les quatre suivantes, renfermees dans le systeme Q, 



28243", 85 = 731939,5z- 13208350z' + 6798,41z"- 1825, 56z'", 

-668486", 70 = -13208350z + 413134432z' - 151992,0z" - 34876,7z" 

245",5870 = 6798,41z-151992,0z' + 71,7381z"-3,7401z // 

42",5434 = -1825,56z-34876,7z'-3,7401z" + 55,0710z" 

ces dernieres equations donnent les suivantes, renfermees dans le systeme dDb, 



_/// 



26833", 55 = 671414,7z- 14364541z' + 6674,43z", 
-695430", = -14364541z + 391046861z'-154360,6z", 
242", 6977 = 6674,43z- 154360,6z' + 71,4841z". 

Enfin j'ai conclu de la les deux equations, renfermees dans le systeme ([E]) : 

4172",95 = 48442z + 48020z, -171455", 2 = 48020z + 57725227z'. 
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5. We now apply this method to an example. For this, I have benefited from the immense work that 
Bouvard has just finished on the movements of Jupiter and Saturn, from which he has constructed extremely 
accurate Tables. He has used all the observations from Bradley and from the astronomers that have followed 
him: he has discussed them again and with the greatest care, which has given him 126 equations for the 
movement of Jupiter in longitude and 129 equations for the movement of Saturn. In these equations, Bouvard 
has introduced the mass of Uranus as unknown. Here are the final equations that he has obtained by la 
methode la plus avantageuse: 



7212".600 = 795938z- 12729398z' + 6788.2z"- 1959.0z'" + 696.13z iV + 2602z v , 

-738297".800 = -12729398z + 424865729z' - 153106.5z" - 39749.1z"' - 5459z"' + 5722z v , 

237".782 = 6788.2z- 153106.5z' + 71.8720z" - 3.2252z"' + 1.2484z"' + 1.337 U v , 

-40".335 = -1959.0z - 39749.1z' - 3.2252z" + 57.19Hz'" + 3.6213z' v + 1.1 128z v , 

-343" .455 = 696. 13z - 5459z' + 1.2484z" + 3.6213z"' + 21.543z' v + 46.310z v , 

- 1002".900 = 2602z + 5722z' + 1 .337 lz" + 1 . 1 128z"' + 46.310z ,v + 129z v . 

In these equations, the mass of Uranus is supposed to be ; the mass of Jupiter is supposed to be 
106709 ' z " * s me P r °duct if the equation of the center by the correction of the periapsis firstly employed 
by Bouvard; z'" is the correction of the equation of the center; z n is the secular correction of the mean 
movement; z 1 is the correction of the epoch of the longitude beginning in 1750. The second of the decimal 
degree is taken for unit. 

With these former equations corresponding to the system (|A|), I have concluded the followings, corre- 
sponding to the system (|B]) : 

27441".68 = 743454z - 128448 14z' + 6761.23z" - 1981.45z"' - 237.97z ,v , 

-693812".58 = - 128448 14z + 4246 11920z' - 153165.81z" - 39798.46z"' -7513.15z ! \ 

248". 1772 = 6761.23z- 153165. 81z' + 71.858U" - 3.2367z"' + 0.7684z ! \ 

-31".6836 = - 1981.45z - 39798.46z' - 3.2367z" + 57. 1815z"' + 3.2218z ! \ 

16".5783 = -237.97z-7513.15z' + 0.7684z" + 3.2218z'" + 4.9181z ,v . 

From these equations, I have drawn the four followings, corresponding to the system Q, 

28243".85 = 731939.5z- 13208350z' + 6798.41z" - 1825.56z'", 

-668486".70 = -13208350z + 413134432z' - 151992.0z" -34876.7z"', 

245".5870 = 6798.4k- 151992.0^ +71. 7381z" -3.740k"', 

42".5434 = -1825.56z-34876.7z'-3.7401z" + 55.07 10z"'; 

these latest equations give us the followings, corresponding to the system ([D]), 

26833".55 = 671414.7z- 14364541z' + 6674.43z", 
-695430".0 = -14364541z + 391046861z'-154360.6z", 
242".6977 = 6674.43z - 154360.6z' + 71.4841z". 

Finally I have concluded from there the two equations, corresponding to the system dE}: 
4172".95 = 48442z + 48020z, -171455" .2 = 48020z + 57725227z'. 
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Je m'arrete a ce systeme, parce qu'il est facile d'en conclude les valeurs du poids P relatives aux deux 
elements z et z' que je desirais particulierement de connaitre. Les formules du n° 3 donnent, pour z, 



P = 



(48020) 2 
48442 - v ' 



57725227 



et, pour z! , 

p _ . „„„„ (48020) 2 



... 57725227- 
2Se'« 2 

Le nombre s des observations est ici 129 et Bouvard a trouve 

Se' {i)2 = 31096; 

on a done, pour z, 



48442 



et, pour z!, 



Les equations precedentes donnent 



log/ 5 = 2,0013595; 
log/ 5 = 5,0778624. 

z ' = -0,00305, 
z = 0,08916. 



La masse de Jupiter est 1()67 09 (1 +z'). En substituant pour z! sa valeur pecedente, cette masse devient 
1()7 q 35 . La masse du Soleil est prise pour unite. La probabilite que l'erreur de z' est comprise dans les limites 
±U est, par le n° 1, 

VP 

l'integrale etant prise depuis u = — U jusqu'a u = U. On trouve ainsi la probabilite que la masse de Jupiter 
est comprise dans les limites 

1 1 1 



du e 



-Pu 2 



1070,35 1001067,09' 

egale a ygggg^y; en sorte qu'il y a un million a tres peu a parier contre un que la valeur 1()7 | ) 35 n'est pas en 
erreur d'un centieme de sa valeur; ou, ce qui revient a fort peu pres au meme, qu'apres un siecle de nouvelles 
observations, ajoutees aux precedentes et dicutees de la meme maniere, le nouveau resultat ne differera pas 
du precedent d'un centieme de sa valeur. 

Newton avait trouve, par les observations de Pound, sur les elongations des satellites de Jupiter, la masse 
de cette planete egale a la 1067 e partie de celle du Soleil, ce qui differe tres peu du resultat de Bouvard. 

La masse d'Uranus est En substituant pour z sa valeur pecedente, cette masse devient La 

probabilite que cette valeur est comprise dans les limites 

1 ± ' 1 



17907 4 19504 

est egale a et la probabilite que cette masse est comprise dans les limites 

1 :± 1 1 



est egale a ^ 



17907 5 19504 

215,6 
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I stop with this system, because it is easy to conclude from it the values of the poids P corresponding to 
the two elements z and z' which I particularly wish to know. The formula from n° 3 give, for z, 



P = 



48442 - 



(48020) 2 " 
57725227 



and, for z', 



P = 



25£ / (')2 



57725227 - 



(48020) : 
48442 



The number s of observations is here 129 and Bouvard has found 

Se' {i)2 = 31096; 



we then have, for z, 
and, for z', 

The former equations give 



logP = 2.0013595; 
log/ 5 = 5.0778624. 

z ' = -0.00305, 
z = 0.08916. 

The mass of Jupiter is 106 7 09 (1 +z'). Replacing z' by its former value, this mass becomes 1()7 | ) 35 . The 
mass of the Sun is taken as unity. The probability that the error in z' is between the limit ±U is, from n° 1, 

^ <due- p » 



VP f 



the integral being taken from u = — U to u = U. We then find that the probability for the mass of Jupiter to 
be between the limits 

1 1 1 



1070.35 1001067.09' 

is equal to {qqqqq',' ; so that there is one million to very few to bet against one that the value WJ l 35 is not 
in error of one hundredth of its value; or, which is more or less the same, that after one century of new 
observations, added to the former and discussed in the same manner, the new result does not differ from the 
former of more than one hundredth of its value. 

Newton had found, from the observations of Pound, on the elongations of Jupiter's satellites, the mass 
of this planet equal to the 1067 th part of the Sun, which differs very few from the result of Bouvard. 

The mass of Uranus is j^^. Replacing for z its former value, this mass becomes The probability 
that this value is between the limits 



1 



1 1 



17907 4 19504 



is equal to and the probability that this mass is between the limits 



1 



■±- 



1 1 



17907 5 19504 



is equal to Hfl . 
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Les perturbations qu'Uranus produit dans le mouvement de Saturne etant peu considerables, on ne doit 
pas encore attendre des observations de ce mouvement une grande precision dans la valeur de sa masse. 
Mais, apres un siecle de nouvelles observations, ajoutees aux precedentes et discutees de la meme maniere, 
la valeur de P augmentera de maniere a donner cette masse avec une grande probabilite que sa valeur sera 
contenue dans d'etroites limites; ce qui sera de beaucoup preferable a l'emploi des elongations des satellites 
d'Uranus, a cause de la difficulte d'observer ces elongations. 

Bouvard, en appliquant la methode precedente aux 126 equations de condition que lui ont donnees les 
observations de Jupiter et en supposant la masse de Saturne egale a 3534 q 8 , a trouve 

z = 0,00620 

et 

log/ 5 = 4, 8856829. 

Ces valeurs donnent la masse de Saturne egale a 35^3, et la probabilite que cette masse est comprise dans 
les limites 

-!-*-!— !_ 

3512,3 100 3534,08 

est egale a yj||| . 

Newton avait trouve par les observations de Pound sur la plus grande elongation du quatrieme satellite 
de Saturne, la masse de cette planete egale a 3^, ce qui surpasse d'un sixieme le resultat precedent. II 
y a des millions de milliards a parier contre un que celui de Newton est en erreur, et Ton n'en sera point 
surpris si Ton considere la difficulte d'observer les plus grandes elongations des satellites de Saturne. La 
facilite d'observer celles des satellites de Jupiter a rendu, comme on l'a vu, beaucoup plus exacte la valeur 
que Newton a conclue des observations de Pound. 
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The perturbations that Uranus induces in the movement of Saturn being negligible, we should not expect 
a great accuracy in the value of the mass from these observations of the movement. But, after a century of 
new observations, added to the previous and discussed in the same manner, the value of P increases so that 
the mass is given with a large probability that its value is contained within tight bounds; which is a lot better 
than using the elongations of the Uranus' satellites, because these elongations are difficult to observe. 

Bouvard, applying the former method to 126 equations given from the observations of Jupiter and as- 
suming that the mass of Saturn is equal to 3534 z 08 , has found 

z = 0.00620 

and 

logP = 4.8856829. 

These values give the mass of Saturn equal to 35^ 3 , and the probability that this mass is between the limits 

-J_±-L_L_ 

3512.3 100 3534.08 

is equal to j]^ • 

Newton has found, from Pound's observations on the largest elongations of Saturn's fourth satellite, that 
the mass of this planet is equal to -J^, which overestimates from than one sixth the former result. There 
are millions of billions to bet against one that the one of Newton is in error, and we should not be surprised 
considering the difficulty to observe the greatest elongations of Saturn's satellites. The easiness to observe 
the ones from Jupiter's satellites has given, as we have seen, a much more exact value than the one concluded 
by Newton from Pound's observations. 
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4 Comments 



Although Laplace presents his algorithm for two variables, in Sections 4.1 and 4.2 we will assume that he 



is only seeking one variable and its standard deviation. This makes explanations easier. Then in Section 4.3 
we generalize to two variables. 



4.1 Laplace's algorithm as a factorization algorithm 

The procedure used by Laplace is a variant of the Cholesky factorization of the normal equations. We do not 
claim that Laplace interprets his algorithm as a factorization. We state that, in Matrix Computation term, we 
can interpret Laplace's algorithm as a factorization. 

In Matlab notation, if we initialize M with the lower part of A r A in input, his algorithm writes 

for k=n:-l:2, 

M(l:k-1, l:k-l) = M ( 1 : k-1 , 1 : k-1 ) - M (1 : k-1, k) *M (1 : k-1, k) ' /M(k,k); 

end 

The operation is a symmertic rank- 1 update and Laplace only updates the lower part of the matrix M at 
each step. After this operation, one obtains a reverse square-root-free Cholesky factorization of the form: 

(A r A) = (D M 1 -M) r M, 

where Dm is the matrix corresponding to the diagonal of M. This is a reverse Cholesky factorization because 
it is a upper triangular matrix times a lower triangular matrix as opposed to being a lower triangular matrix 
times a upper triangular matrix. It is square root free because the left factor, (D^ 1 • M) , has a unit diagonal 

1/2 

and the right factor, M, has a non-unit diagonal, reverse Cholesky (with square roots) gives to each 
factor so that the factorization writes 

(A^A) = (D M 1 / 2 .M) r (D M 1 / 2 .M)=L^L. 



Cholesky factorization reverse square-root-free Cholesky factorization 




In Matlab notation, after Laplace's algorithm, we could compute the solution with a backward and a 
forward solve. 

z = (diag(diag(M) ) \M) ' \ alpha; 
z = M \ z; 

where the first line is the backward solve and the second line is the forward solve. 
In Laplace's algorithm, the backward solve is done on the fly as 

z = alpha; 
for k=n:-l:2, 

z(l:k-l) = z(l:k-l) - M (1 : k-1, k) *z (k) / M(k,k); 

end 
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On the fly means that this loop is inserted in the loop of the factorization. Since Laplace only wishes the 
first variable, the forward solve reduces to 



z(l) = z(l) / M(l,l); 

Another way to understand Laplace's factorization is to follow Trefethen and Bau's explanations of 
LU 01. Laplace is applying a sequence of unit upper triangular matrices, U,, to A r A in order to reduce 
A r A to lower triangular form. We obtain 



/i 



V 



U n _i 



i * 
i 



i * 
i 



i 



u 2 



V 



/ 



V 



And with a few "strokes of luck", we can prove that 



A^A 




(U J ,_i,...,U 2 ,Ui 



(Dm ! -M) j 



so that we finally recover our factorization: (A r A) = (D^ 1 • M) M. 

If we follow this interpretation then, the backward solve is the application of the sequence of unit upper 
triangular matrices to the right-hand side itself. 

A third way to understand Laplace's factorization is to follow Laplace's explanations. Laplace is im- 
plicitly computing the L factor of the QL factorization of A from the normal equations A T A. 



4.2 Laplace's algorithm to compute the standard deviation of a variable 



We have seen in Section [2.2.2| fhat we can derive quite easily the standard deviation of the first variable of a 
statistical model from the QL factorization of the regression matrix. The relation is given in Equation Q. 

Since the L factor obtained by the QL factorization is also the reverse Cholesky factor, Equation ([5]) 
explains how to compute the standard deviation of the first variable with reverse Cholesky. We think this is 
the best way to understand Laplace's algorithm with our contemporary tools. 

4.3 Laplace's algorithm to compute the standard deviation of two variables 

Laplace indeed is interested in the standard deviation of the two first variables. In this case, he stops his 
reduction process at the 2-by-2 matrix corresponding to 

A r A =( P«-2P«-2 Pn- 2 qn-2 \ = ( 48442 48020 
n - 2 n ~ 2 V €-iVn-2 <£_ 2 qn-2 J V 48020 57725227 

The covariance matrix of u and u' is then given by 

t a .N-i _ p£_2<l»-2 V ! _ 48442 48020 N ' 



a a (A„_ 2 A„_ 2 ) _ a^ q ,_ 2pn 2 ^ J - ° b y 48020 57725227 

2 1 / 57725227 -48020 

= of 1 



fc 48442 * 57725227 - 48020 2 V -48020 48442 
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1 1 1 / 1 1 2 o 

If we use the fact that -j— ||e || is an unbiased estimate of <5 b and approximate the term s — n with s, we 
get that the standard deviation of z, the first variable, and the standard deviation of z', the second variable, 
are equal to 



1 



57725227 



48442 * 57725227 - 48020 2 



and c' 2 



1 



48442 



48442 * 57725227 - 48020 2 ' 



Laplace gives these two formulae for the poids, P, and not for the standard deviation. Using the relation 
P = (2g 2 ) -1 we find the equation of Laplace for z (top of p 24 of this document) 



2Se'(')2 



48442 - 



(48020) 2 " 
57725227 



and, for z', 



2S£ / « 2 



57725227 - 



(48020) : 
48442 



5 Numerical example 

There are two problem sets for Laplace to apply his algorithm. The first one computes the mass of Jupiter, the 
second one computes the mass of Saturn. Laplace uses observations from Bouvard where (s = 126, n = 6) 
for Jupiter and (s = 129, n = 6) for Saturne. Actually Laplace only uses from Bouvard the 6-by-6 normal 
equations and the norm of the residual of the least squares problem, e'. On the 6 unknowns (in the x vector), 
Laplace only seeks one, the second variable z' ■ The mass of Jupiter in term of the mass of the Sun is given 
by z' and the formula: 

mass of Jupiter = . 

F 1067.09 

It turns out that the first variable, z, represents the mass of Uranus through the formula 

mass of Uranus = — . 

19504 

Same approach holds for Saturn, so Laplace will indeed compute and report the mass of three planets in his 
manuscript. 

Note that at this time, Bouvard knew that he did not understand the behavior of Uranus. He conjectured 
that another planet should exist to explain the anomality in the observed behavior of Uranus. The mass of 
Uranus is introduced as the auxiliary variable z to try to cure the problem. Laplace correctly predicts that the 
computed mass for Uranus is not reliable. For the anecdoct, the missing planet was Neptun and was found 
by Johann Gottfried Galle three years after the death of Bouvard. 

The number of operations performed by Bouvard is quite remarkable. For the computation of the mass of 
Jupiter, Laplace accredited Bouvard for the computation of the normal equations (A r A) and of the residual 
norm (e ; = Ax — b), this makes about sn 2 + 2sn operations. For this numerical example, Laplace performed 
the Cholesky factorization which is about n 3 /3. This represents 6,048 operations for Bouvard and a mere 
72 for Laplace! For the computation of the mass of Saturn, the comparison is even worse since Bouvard 
performed all the operations and reported the results to Laplace. We note that this means that Laplace has 
explained his algorithm to Bouvart. 

The computation of Laplace proved to be quite exact. In Table [T] we compare them with the current 
NASA values. We see that the values for Jupiter and Saturn in Laplace are quite close from the NASA 
ones. The value for Uranus is quite far as can be expected from its large variance. (Laplace would have said 
its small poids.) We also note that the NASA values are within Laplace's bounds for Saturn and Uranus. 
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Planet 


Laplace/Bouvard 


estimation of the probability of error 


NASA 


Jupiter 


(1059) 1070 (1081) 


1/1,000,000 


1048 


Uranus 


(14564) 17918 (23241) 


1/2,509 


22992 


Saturn 


(3477) 3512 (3547) 


1/11,328 


3497 



Table 1: Fraction of the mass of the Sun. The computed values from Bouvard are given in bold and the 
bound from Laplace in parenthesis. Laplace proved that his value for the mass of Uranus was not reliable 
(He was right.) The interval of confidence for Uranus and Saturn from Laplace are correct (i.e. the NASA 
values are in these intervals). 



The value for Jupiter is not within Laplace's bound which means that the noise in the observations was not 
normal. 

In Table [2j we perform the computation of Laplace again using 64-bit arithmetic and we report the 
incorrect digits in his computation. It is interesting to see that Laplace conserves a fix number of significant 
digits along the computation. We can therefore say that Laplace was computing in floating-point arithmetic. 

We note that, while the condition number of A r A is failry large (above 10 8 ), we can equilibrate the 
matrix A r A with a diagonal scaling S equal to the inverse of the square-root of the diagonal elements. In 
this case, the scaled normal equations matrix, S(A T A)S, has ones on its diagonal and its condition number 
of 104. So, up to a diagonal scaling, the system that Laplace is considered is well-conditioned. 

We can check the value given by Laplace for the variance of the variable z and z' ■ On the one hand, 
Laplace gives the poids of z as log/ 3 = 2.0013595 so we obtain that the standard deviation of z is given 
by 1/sqrt (2) /sqrt (10" (2 .0013595) ) that is a z = 0.0706. On the other hand we can use the standard 

formula a z = c fo gentry ( 1, 1) of (A r A) _1 . In Matlab, this gives sqrt (31096/129) *sqrt ( invATA (1,1) ) 
we obtain a z = 0.0707. For the variable z' Laplace gives its poids as 5.0778624, therefore o' z = 0.002044343. 
Directly from the normal equations, we would have found o z = 0.002044348. 

Laplace intreprets his result by giving an interval with a confidence level. For example, once, Laplace 
has computed z' = 0.08916, the variable such that the mass of Jupiter is l ^ J z 09 , and its associted poids 
(P = 105.0778624^ Lap i ace uses t h e fact that 

1/100 , -pJ ioooooo 

du e w 

1000001 

to claim that there is one chance out of one million for the computed value of z' to be between — 1/100 and 
1/100 of its exact value. This means that there is one chance out of one million for the mass of Jupiter to be 
between 1/1,081 and 1/1,059 the one of the Sun. 

In the same manner, once, Laplace has computed z = —0.00305, the variable such that the mass of 

l+z 

17907 ' 




Uranus is . kfe , and its associted poids (P = 1Q 2 0013595 ), Laplace uses the fact that 




_ Ph 2_ 2508 
due « — — - 
2509 

to claim that there is one chance out of 2, 509 for the computed value z to be between — 1/4 and 1/4 of its 
exact value. This means that there is one chance out of 2, 509 for the mass of Uranus to be between 1/23,241 
and 1/14,564 the one of the Sun. 
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STEP B: 


743454 -12844814 6761.23 -1981.45 -237.97 
424611920 -153165.81 -39798.46 -7513.15 
71.8581 -3.2367 0.7684 
57.1815 3.2218 
4.918 


27441.68 
-693812.58 
248.1772 
-31.6836 
16.5783 


743454 -12844814 6761.23 -1981.45 -237.97 
424611920 -153165.81 -39798.46 -7513.15 
71.8581 -3.2367 0.7684 
57.1815 3.2218 
4.918 


27441.64 
-693812.58 
248.1772 
-31.6836 
16.5783 



STEP A: 

795938 -12729398 
424865729 



6788.2 
-153106.5 
71.8720 



-1959.0 
-39749.1 
-3.2252 
57.1911 



696.13 
-5459 
1.2484 
3.6213 
21.543 



2602 
5722 
1.3371 
1.1128 
46.310 
129 



7212.600 
-738297.800 
237.782 
-40.335 
-343.455 
-1002.900 



STEP C: 



731939.5 



-13208350 
413134432 



6798.41 
-151992.0 
71.7381 



-1825.56 
-34876.7 
-3.7401 
55.0710 



28243.85 
-668486.70 
245.5870 
42.5434 



731939.2 



-13208360 
413134201 



6798.41 
-151991.9 
71.7380 



-1825.55 
-34876.6 
-3.7401 
55.0709 



28243.86 
-668486.18 
245.5870 
42.5441 



STEP D: 



671414.7 



-14364541 
391046861 



6674.43 
-154360.6 
71.4841 



26833.55 
-695430.0 
242.6977 



671423.6 



-14364485 
391046869 



6674.43 
-154360.6 
71.4841 



26833.57 
-695429.6 
242.6977 



STEPE: 



48442 48020 
57725227 

48227 48021 
57725258 



4172.95 
-171455.2 

4173.00 
-171355.9 



STEP F: 



Zo = -0.00305 
zi =0.08916 



zo = -0.00304 
zi =0.08916 



Table 2: Laplace computations. Comparison with "exact" computation, ("exact" means 64-bit arithmetic.) 
First line is Laplace's value. Second line is the value computed with 64-bit arithmetic and Laplace's algo- 
rithm. 
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